This invention relates to systems and methods for reconstructing acoustic and/or electromagnetic properties of two or three-dimensional objects using diffraction tomographic procedures. The systems and methods are useful in the medical diagnostics arts such as ultrasonic and X-ray tomography, as well as in the geophysical prospecting arts, such as well-to-well tomography and subsurface seismic prospecting.
In conventional parallel beam transmission computed tomography, an object profile O(x,y) representing, for example, an X-ray attenuation coefficient, is reconstructed from the object's projections EQU P.sub..theta. (.xi.)=.intg.O(x,y)d.eta., (1)
where (.xi.,.eta.) denote the coordinates in a cartesian coordinate system rotated by the angle .theta. relative to the (x,y) system as shown in FIG. 1a. The reconstruction of O(x,y) from P.sub..theta. (.xi.) (which is called the "reconstruction of the object" or the "reconstruction of the object profile" by those skilled in the art) is made possible by the projection-slice theorem which states that the one-dimensional Fourier transform of P.sub..theta. (.xi.) is a slice at angle .theta. through the two-dimensional Fourier transform O(K.sub.x,K.sub.y) of O(x,y) as shown in FIG. 1b. Thus, the projections P.sub..theta. (.xi.) yield, upon Fourier transformation, an ensemble of slices of the two-dimensional Fourier transform O(K.sub.x,K.sub.y) of O(x,y). O(x,y) can then readily be reconstructed from this ensemble of slices via its Fourier integral representation expressed in circular polar coordinates. In particular, one finds that ##EQU1## where ##EQU2## is the one-dimensional Fourier transform of P.sub..theta. (.xi.). In Equation (2), the K integration is limited to the interval -W to W since, in any application, P.sub..theta. (K) will be known only over a finite bandwidth which is taken to be -W to W. Because of this restriction, the resulting reconstruction is a low pass filtered version of O(x,y).
In practice, in place of the Fourier integral representation given in Equation (2), a technique known in the art as the filtered backprojection algorithm is used to reconstruct O(x,y). Defining the filter function ##EQU3## which, when convolved with a low pass filter, is known in the art as the standard X-ray tomography deblurring filter, it can be shown (see, for example, A. C. Kak, "Computerized Tomography with X-ray Emission and Ultrasound Sources"; Proc. IEEE 67, pp. 1245-1272 (1974)) that Equation (2) can be written in the form ##EQU4## where Q.sub..theta. (t), called "filtered projections", are given by ##EQU5## The process of reconstructing O(x,y) according to the prior art thus consists of first filtering the projections with the filter function h(t) and then backprojecting the filtered projections.
The process of backprojection used in the prior art consists of assigning any given pixel value (x,y) within the image array the value of the filtered projection Q.sub..theta. (t) at the point EQU t=x cos .theta.+y sin .theta.. (7)
Backprojection at any given view angle .theta. thus results in a partial image consisting of parallel straight lines as defined by Equation (7) with the grey level of each line being assigned the value Q.sub..theta. (t). Following the process of backprojection at a fixed view angle, the entire process (filtering followed by backprojection) is repeated at different angles and the resulting partial image arrays are added pixel by pixel onto the previously existing partial image (hence, the integral over d.sub..theta. in Equation 5). In this manner, the profile of the object is reconstructed. In practice, of course, the filtered projection is only known for a finite number of angles so that a discrete approximation (e.g. a summation) to the backprojection integral Equation (5) is employed.
The success of transmission computed tomography clearly depends on the availability of the projections P.sub..theta. (.xi.). In X-ray tomography, these projections are available in the form of measurements of the log amplitude [1/2 log (intensity)] of the electromagnetic field produced when a plane, monochromatic X-ray beam propagating along the .eta. axis is incident to the object as is illustrated in FIG. 2a. At X-ray wavelengths (.apprxeq.1 .ANG.) the incident plane wave suffers very little scattering (refraction or diffraction) in passing through the object. Absorption does occur, however, so that the negative of the log amplitude of the field along the line .eta.=1.sub.0 is simply the integrated value of the X-ray absorbtion coefficient .mu.(x,y) along the travel path (i.e. along the straight line .xi.=constant). The object profile O(x,y) in this case is then .mu.(x,y) and the projections P.sub..theta. (.xi.) are simply the negative of the log amplitude of the electromagnetic field as measured along the line .eta.=1.sub.0.
In ultrasonic tomography, the wavelength is considerably greater (.apprxeq.1 mm.) than that of X-rays so that an acoustic wave experiences scattering, in the form of refraction and diffraction, in the process of propagating through an object. The log amplitude of the acoustic field along the line .eta.=1.sub.0 resulting from an insonifying plane wave propagating along the .eta. axis is then not simply the projection of the acoustic absorbtion coefficient along the line .xi.=constant as is the case in X-ray tomography. This breakdown in the relationship between measured log amplitude and projections of the absorbtion coefficient is one of the major reasons why conventional computed tomography is of limited use in ultrasonic applications.
Although it is not possible, in general, to obtain an exact expression for the acoustic field generated by the interaction of an insonifying plane wave with an acoustic (diffracting) object O(x,y), expressions within the so-called first Born and Rytov approximations are available.
The Born and Rytov approximations are essentially "weak scattering" approximations and, as such, will be valid for cases where the deviations of the object profile O(x,y) from zero are small. Which approximation, Born or Rytov, is employed depends upon the specific application, although the Rytov approximation is generally thought to be superior in medical ultrasound tomography (See, for example, M. Kaveh, M. Soumekh and R. K. Mueller, "A Comparison of Born and Rytov Approximations in Acoustic Tomography", Acoustical Imaging, Vol. 11 Plenum Press New York, (1981).) A more complete treatment of the Born and Rytov approximations is found in A. J. Devaney, "Inverse Scattering Theory Within the Rytov Approximation," Optics Letters 6, pp. 374-376 (1981). Within these approximations, it has been shown by the inventor in "A New Approach to Emission and Transmission CT" 1980 IEEE Ultrasonics Symposium Proceedings, pp. 979-983 (1980); "Inverse Scattering Theory Within the Rytov Approximation" Opics Letters 6, pp. 374-376 (1981); and by R. K. Mueller, et al., "Reconstructive Tomography and Application to Ultrasonics", Proc. IEEE 67, pp. 567-587 (1979), that for an incident monochromatic plane wave propagating along the .eta. axis, the profile O(x,y) is related to the processed signals D.sub..theta..sbsb.0 (.xi.,.omega.) (hereinafter referred to as "data") taken along the line .eta.=1.sub.0 via the equation ##EQU6## where the variable .kappa. can assume all values in the range[-k,k]. Here .theta..sub.0 =.theta.+.pi./2 is the angle that the wave vector of the insonifying plane wave makes with the x axis, .omega. is the angular frequency, and k=2.pi./.lambda. is the wavenumber of the insonifying acoustic field. The data D.sub..theta..sbsb.0 (.xi.,.omega.), are shown in the case of the Born approximation, to be simply related to the complex amplitude of the scattered field portion of the acoustic field evaluated along the line .eta.=1.sub.0 and, in the case of the Rytov approximation, to the difference in the complex phase of the total field and incident field evaluated along this same line.
For discussion regarding diffraction tomography, it is convenient to introduce two unit vectors; ##EQU7## where .xi. and .eta. are unit vectors along the .xi., .eta. axes. The unit vector s.sub.0 is simply the unit propagation vector of the insonifying plane wave. For values of .vertline..kappa..vertline.&lt;k, Equation (8) can be expressed in terms of s and s.sub.0 in the form ##EQU8## From Equation (10) it may be concluded that the one-dimensional Fourier transform D.sub..theta..sbsb.0 (.kappa.,.omega.) of D.sub..theta..sbsb.0 (.xi.,.omega.) is equal to the two-dimensional Fourier transform O(K.sub.x,K.sub.y) of the profile O(x,y) evaluated over the locus of points defined by EQU K.ident.K.sub.x X+K.sub.y y=k(s-s.sub.0), (11)
where the unit vector s assumes all values for which s.multidot.s.sub.0 &gt;0. Equation 10 may thus be seen as the diffraction tomography equivalent of the projection slice theorem. The K values satisfying Equation (11) for fixed s.sub.0 with s.multidot.s.sub.0 &gt;0, are seen to lie on a semicircle, centered at -ks.sub.0 and having a radius equal to k. By changing the direction of propagation of the insonifying plane wave, it is thus possible to sample O(K.sub.x,K.sub.y) over an ensemble of semicircular arcs such as illustrated in FIG. 2b. The ultrasonic tomographic reconstruction problem then consists of estimating O(x,y) from specification of O(K.sub.x,K.sub.y) over this ensemble of arcs.
The problem of reconstructing an object profile from specification of its Fourier transform over an ensemble of circular arcs such as shown in FIG. 2b is an old one, with proposed solutions having first occurred in X-ray crystallography, H. Lipson and W. Cochran, The Determination of Crystal Structures (Cornell University Press, Ithica, N.Y., 1966), and later in inverse scattering applications using both X-rays, B. K. Vainshtein, Diffraction of X-rays by Chain Molecules (John Wiley, New York, 1977), and visible light, E. Wolf, "Three-dimensional Structure Determination of Semi-Transported Objects from Holographic Data", Optics Comm., pp. 153-156 (1969). More recently, these inverse scattering solutions have arisen in ultrasonics, (see e.g. S. J. Norton and M. Linzer, "Ultrasonic Reflectivity Imaging in Three Dimensions: Exact Inverse Scattering Solutions for Phase, Cylindrical and Spherical Apertures," IEEE Trans Biomed. Eng. BME-28, pp. 202-220 (1981)) and also in ultrasonic tomography which has become known as diffraction tomography to distinguish it from conventional tomography where diffraction effects are not taken into account. In all of these applications, the reconstruction of the profile O(x,y) is made possible by the fact that by varying s.sub.0 over the unit circle, one can determine O(K.sub.x,K.sub.y) at any point lying within a circle, centered at the origin and having a radius equal to .sqroot.2k. A low pass filtered version O.sub.LP (x,y) of O(x,y) can then be obtained by means of the Fourier integral representation ##EQU9## where vector notation is used such that r=(x,y), K=(K.sub.x,K.sub.y) and d.sup.2 K stands for the differential surface element in K.
In any application, the profile O(r) will vanish beyond some finite range so that a Fourier series expansion for O.sub.LP (r) can be used in place of the Fourier integral (12). According to the solutions proposed in the art, in the Fourier series expansion, the transform O(K) must then be known for K values lying on a regular, square sampling grid in K space. To obtain the requisite sample values therefore requires that either one interpolate O(K.sub.x,K.sub.y) from the sample values given on arcs over which this quantity is known (interpolation method) or alternatively, that one select propagation directions s.sub.0 such that at least one arc falls on every required sample point in the array (direct sampling method). Discussion and examples of the interpolation method can be found in W. H. Carter, "Computerized Reconstruction of Scattering Objects from Holograms", J. Opt. Soc Am 60., pp. 306-314 (1976), while the direct sampling method procedure is followed in e.g. A. F. Fercher et al., "Image Formation by Inversion of Scattered Field Data: Experiments and Computational Simulation", Applied Optics 18, pp. 2427-2431 (1979).
Both the interpolation and direct sampling methods for obtaining sample values of O(K) over a regular square sampling grid have serious drawbacks. Interpolation of O(K) from sample values given on semicircular arcs suffers from the same limitations as does interpolation from slices in K space. Reconstructions based on two-dimensional Fourier Transforms from interpolated sample values in conventional tomography are well known to be inferior in accuracy to those obtained via filtered backprojection. Direct sampling methods, on the other hand, are undesirable because they generally require an extremely large number of s.sub.0 values (i.e., many experiments). In addition, both procedures place severe if not insurmountable demands on the experimentalist in that they require extreme accuracy in registering data collected from repeated tomographic procedures.
The present invention is a great improvement over the prior art in that a method and system are provided which allow the low pass filtered object profile O.sub.LP (r) to be reconstructed directly from processed signals D.sub..theta..sbsb.0 (.xi.,.omega.) without the necessity of first determining O(K) over a regular square sampling grid in K space as was required in the prior art. The provided methods and systems employ the so-called filtered backpropagation algorithm (so named by the inventor) which is closely analogous to the filtered backprojection algorithm of conventional transmission tomography. The two techniques differ in one important respect: in the backprojection process the filtered projections Q.sub..theta. (t) are continued back through the object space along parallel straight line paths; whereas in the process of backpropagation, the filtered phases are continued back in the manner in which they were diffracted. In other words, in the backprojection process, the points along a wave of energy may each be considered to follow straight parallel paths through the object. In the backpropagation process, the points are each considered to undergo scattering throughout the object. For example; the Rytov and Born approximations permit an analysis of the problem by suggesting that only first scatterings be taken into account, and that secondary scatterings may be ignored as insignificant. Thus, in backpropagation, quantities Q(t), related to the data D.sub..theta..sbsb.0 (.xi.,.omega.) through the filtering operation of equation (6), are mapped onto the object space via an integral transform that is the inverse of the transform that governs phase propagation within the Rytov approximation. It is natural then to call this process backpropagation since it corresponds mathematically to the inverse of the usual forward propagation process. The relationship which is employed and embodied in the filtered backpropagation technique has been conveniently denoted by the inventor as the filtered backpropagation algorithm to emphasize its formal similarity to the filtered backprojection algorithm of conventional tomography.
Thus, for purposes of this application and the claims herein, the term "backpropagation" shall be defined to mean that operation that is the inverse or approximate inverse of a forward propagation process. The term "filtered backpropagation technique" shall be defined to describe any diffraction tomographic technique for the partial or complete reconstruction of an object where a filtered real or complex amplitude and/or filtered real or complex phase of a wave is backpropagated into the object space; i.e., is propagated back into object space according to the inverse or approximate inverse of the way in which the wave was originally diffracted. The filtered backpropagation technique is usually implemented in the form of a convolution of filters. For purposes of brevity, such an implementation will identically be called the filtered backpropagation technique. A "backpropagation filter" shall be defined to describe that filter of the filtered backpropagation technique which accounts for diffraction in the backpropagation of the phase; e.g. in ultrasound diffraction tomography, the filter of the filtered backpropagation technique which is not the standard X-ray tomographic filter (or a variation thereon). A "filtered backpropagation operation" shall be defined as any procedure which employs the filtered backpropagation technique. Also, for purposes of this application and the claims herein, it should be understood that the waves of energy which propagate and diffract according to the invention include but are not limited to sonic or electromagnetic waves. The term "sonic wave" shall be interpreted as broadly as possible and shall be understood to include all elastic wave phenomena in liquid and solid materials including, but not limited to, acoustic, compressional, shear, and elastic waves. The term "acoustic wave" shall be interpreted herein to be the equivalent of "sonic wave". The term "electromagnetic wave" shall also be interpreted in its broadest sense and shall include, but not be limited to infrared rays, X-rays, and the class known as "optics".